library("rstudioapi")     
##Set working directory
setwd(dirname(getActiveDocumentContext()$path))
getwd()
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(estimatr)
library(tidyverse)
library(geosphere)
library(DescTools)
require(dplyr)

load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

# we need 6 version of each model: 
working_final = final_finalV2
working_finalA0 = working_final[which(working_final$Attended == 0),]
working_finalA1 = working_final[which(working_final$Attended == 1),]
#working_finalF = working_final[which(working_final$Q2.3 == "Female"),]
#working_finalM = working_final[which(working_final$Q2.3 == "Male"),]
working_finalW0 = working_final[which(working_final$WhatsApp == 0),]
working_finalW1 = working_final[which(working_final$WhatsApp == 1),]


### Model 7c combinations: 
L_7cG <- glm.cluster(vaccine_intention ~ cash_cat + Gender + Age + 
                     Income + Employed + Education + cash_cat*Gender, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_7cA <- glm.cluster(vaccine_intention ~ cash_cat + Gender + Age + 
                     Income + Employed + Education + cash_cat*Education, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_7cW <- glm.cluster(vaccine_intention ~ cash_cat + Gender + Age + 
                     Income + Employed + Education + cash_cat*WhatsApp, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)

tab_7cA = exp(cbind("Odds ratio" = coef(L_7cA), confint.default(L_7cA, level = 0.95)))
rownames(tab_7cA)[1:10] = c("Intercept", "Cash", "Health", "Male", "Age", "Mean-Food", "Employed", "High-Educate", "Low-Educate", "Medium-Educate")
tab_7cW = exp(cbind("Odds ratio" = coef(L_7cW), confint.default(L_7cW, level = 0.95)))
rownames(tab_7cW)[1:10] = c("Intercept", "Cash", "Health", "Male", "Age", "Mean-Food", "Employed", "High-Educate", "Low-Educate", "Medium-Educate")
tab_7cG = exp(cbind("Odds ratio" = coef(L_7cG), confint.default(L_7cG, level = 0.95)))
rownames(tab_7cG)[1:10] = c("Intercept", "Cash", "Health", "Male", "Age", "Mean-Food", "Employed", "High-Educate", "Low-Educate", "Medium-Educate")

col_7cA = c()
for (i in c(1:dim(tab_7cA)[1])){
  col_7cA[i] = as.character(paste(format(round(tab_7cA[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_7cA[i,2],2), nsmall = 2), format(round(tab_7cA[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_7cA = data.frame(col_7cA);rownames(col_7cA) = rownames(tab_7cA);colnames(col_7cA) = "ModelA"

col_7cW = c()
for (i in c(1:dim(tab_7cW)[1])){
  col_7cW[i] = as.character(paste(format(round(tab_7cW[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_7cW[i,2],2), nsmall = 2), format(round(tab_7cW[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_7cW = data.frame(col_7cW);rownames(col_7cW) = rownames(tab_7cW);colnames(col_7cW) = "ModelW"

col_7cG = c()
for (i in c(1:dim(tab_7cG)[1])){
  col_7cG[i] = as.character(paste(format(round(tab_7cG[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_7cG[i,2],2), nsmall = 2), format(round(tab_7cG[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_7cG = data.frame(col_7cG);rownames(col_7cG) = rownames(tab_7cG);colnames(col_7cG) = "ModelG"


cases7c = merge(col_7cA, col_7cW, by = 0, all = TRUE)
rownames(cases7c) = cases7c[,1]
cases7c[,1] = NULL
cases7c = merge(cases7c, col_7cG, by = 0, all = TRUE)
rownames(cases7c) = cases7c[,1]
cases7c[,1] = NULL

output = cases7c[c(2,14,1,18,15,20,17,13,19,21,16,3,5,4,8,10,9,6,11,7,12),]
#stargazer(output, type = "text", summary = FALSE, out="Tables/Table17.txt")
#stargazer(output, type = "latex", summary = FALSE, out="Tables/Table17.tex")

no.obs = c((L_7cA$glm_res$df.null), (L_7cW$glm_res$df.null), (L_7cG$glm_res$df.null))
aic_mod = round(c(L_7cA$glm_res$aic, L_7cW$glm_res$aic, L_7cG$glm_res$aic),2)
log_lik = round(c(as.vector(logLik(L_7cA$glm_res)), as.vector(logLik(L_7cW$glm_res)), as.vector(logLik(L_7cG$glm_res))),2)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(output)

# merge output with mod spec: 
output2 = rbind(output, mod_spec)
stargazer(output2, type = "text", summary = FALSE, out="ExtData_Table6_Part1.txt")
stargazer(output2, type = "latex", summary = FALSE, out="ExtData_Table6_Part1.tex")

part1 = output2

# Model 8c: 
L_8cG <- glm.cluster(vaccine_reported_combo ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*Gender, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_8cA <- glm.cluster(vaccine_reported_combo ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*Education, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_8cW <- glm.cluster(vaccine_reported_combo ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*WhatsApp, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)

tab_8cA = exp(cbind("Odds ratio" = coef(L_8cA), confint.default(L_8cA, level = 0.95)))
rownames(tab_8cA)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")
tab_8cW = exp(cbind("Odds ratio" = coef(L_8cW), confint.default(L_8cW, level = 0.95)))
rownames(tab_8cW)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")
tab_8cG = exp(cbind("Odds ratio" = coef(L_8cG), confint.default(L_8cG, level = 0.95)))
rownames(tab_8cG)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")

col_8cA = c()
for (i in c(1:dim(tab_8cA)[1])){
  col_8cA[i] = as.character(paste(format(round(tab_8cA[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_8cA[i,2],2), nsmall = 2), format(round(tab_8cA[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_8cA = data.frame(col_8cA);rownames(col_8cA) = rownames(tab_8cA);colnames(col_8cA) = "ModelA"

col_8cW = c()
for (i in c(1:dim(tab_8cW)[1])){
  col_8cW[i] = as.character(paste(format(round(tab_8cW[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_8cW[i,2],2), nsmall = 2), format(round(tab_8cW[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_8cW = data.frame(col_8cW);rownames(col_8cW) = rownames(tab_8cW);colnames(col_8cW) = "ModelW"

col_8cG = c()
for (i in c(1:dim(tab_8cG)[1])){
  col_8cG[i] = as.character(paste(format(round(tab_8cG[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_8cG[i,2],2), nsmall = 2), format(round(tab_8cG[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_8cG = data.frame(col_8cG);
rownames(col_8cG) = rownames(tab_8cG)
colnames(col_8cG) = "ModelG"

cases8c = merge(col_8cA, col_8cW, by = 0, all = TRUE)
rownames(cases8c) = cases8c[,1]
cases8c[,1] = NULL
cases8c = merge(cases8c, col_8cG, by = 0, all = TRUE)
rownames(cases8c) = cases8c[,1]
cases8c[,1] = NULL

output = cases8c[c(2,14,1,18,15,20,17,13,19,21,22,16,3,5,4,8,10,9,6,11,7,12),]
#stargazer(output, type = "text", summary = FALSE, out="Tables/Table20.txt")
#stargazer(output, type = "latex", summary = FALSE, out="Tables/Table21.tex")

no.obs = c((L_8cA$glm_res$df.null), (L_8cW$glm_res$df.null), (L_8cG$glm_res$df.null))
aic_mod = round(c(L_8cA$glm_res$aic, L_8cW$glm_res$aic, L_8cG$glm_res$aic),2)
log_lik = round(c(as.vector(logLik(L_8cA$glm_res)), as.vector(logLik(L_8cW$glm_res)), as.vector(logLik(L_8cG$glm_res))),2)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(output)

# merge output with mod spec: 
output2 = rbind(output, mod_spec)
stargazer(output2, type = "text", summary = FALSE, out="ExtData_Table6_Part2.txt")
stargazer(output2, type = "latex", summary = FALSE, out="ExtData_Table6_Part2.tex")

part2 = output2

rownames(part2)
rownames(part1)

tb1 = merge(part1, part2, by = 0, all = TRUE);rownames(tb1) = tb1[,1];tb1[,1] = NULL

tb2 = tb1[c("Cash","Health","Age" ,"Male","High-Educate","Medium-Educate","Low-Educate","Employed","Mean-Food","Social Media", "WhatsApp", "Intercept",     
      "cash_catcash:EducationHigh","cash_catcash:EducationMedium", "cash_catcash:EducationLow","cash_catCDC:EducationHigh",   
      "cash_catCDC:EducationMedium", "cash_catCDC:EducationLow", "cash_catcash:GenderMale","cash_catCDC:GenderMale",      
      "cash_catcash:WhatsApp", "cash_catCDC:WhatsApp", "Observations", "Akaike Inf. Crit.", "Log Likelihood")
    ,]

#stargazer(tb2, type = "text", summary = FALSE, out="ExtData_Table6_FULL.txt")
stargazer(tb2, type = "latex", summary = FALSE, out="ExtData_Table6_FULL.tex")






# Model 9c: 
L_9cG <- glm.cluster(ActVacApril ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*Gender, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_9cA <- glm.cluster(ActVacApril ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*Education, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)
L_9cW <- glm.cluster(ActVacApril ~ cash_cat  + Gender + Age + 
                     Education + Income + Employed + SNMetric + cash_cat*WhatsApp, 
                     data = working_final, family=binomial,
                     cluster = working_final$Q123)

tab_9cA = exp(cbind("Odds ratio" = coef(L_9cA), confint.default(L_9cA, level = 0.95)))
rownames(tab_9cA)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")
tab_9cW = exp(cbind("Odds ratio" = coef(L_9cW), confint.default(L_9cW, level = 0.95)))
rownames(tab_9cW)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")
tab_9cG = exp(cbind("Odds ratio" = coef(L_9cG), confint.default(L_9cG, level = 0.95)))
rownames(tab_9cG)[1:11] = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media")

col_9cA = c()
for (i in c(1:dim(tab_9cA)[1])){
  col_9cA[i] = as.character(paste(format(round(tab_9cA[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_9cA[i,2],2), nsmall = 2), format(round(tab_9cA[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_9cA = data.frame(col_9cA);rownames(col_9cA) = rownames(tab_9cA);colnames(col_9cA) = "ModelA"

col_9cW = c()
for (i in c(1:dim(tab_9cW)[1])){
  col_9cW[i] = as.character(paste(format(round(tab_9cW[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_9cW[i,2],2), nsmall = 2), format(round(tab_9cW[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_9cW = data.frame(col_9cW);rownames(col_9cW) = rownames(tab_9cW);colnames(col_9cW) = "ModelW"

col_9cG = c()
for (i in c(1:dim(tab_9cG)[1])){
  col_9cG[i] = as.character(paste(format(round(tab_9cG[i,1],2), nsmall = 2),paste("(",paste(format(round(tab_9cG[i,2],2), nsmall = 2), format(round(tab_9cG[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
}
col_9cG = data.frame(col_9cG);
rownames(col_9cG) = rownames(tab_9cG)
colnames(col_9cG) = "ModelG"

cases9c = merge(col_9cA, col_9cW, by = 0, all = TRUE)
rownames(cases9c) = cases9c[,1]
cases9c[,1] = NULL
cases9c = merge(cases9c, col_9cG, by = 0, all = TRUE)
rownames(cases9c) = cases9c[,1]
cases9c[,1] = NULL

output = cases9c[c(2,14,1,18,15,20,17,13,19,21,22,16,3,5,4,8,10,9,6,11,7,12),]
#stargazer(output, type = "text", summary = FALSE, out="ExtData_Table6.txt")
#stargazer(output, type = "latex", summary = FALSE, out="ExtData_Table6.tex")
